Quantum algorithm to calculate electromagnetic scattering cross sections 
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We describe a quantum algorithm that generalizes the quantum linear system algorithm [Harrow 
et at, Phys. Rev. Lett. 103, 150502 (2009)] to arbitrary problem specifications. We develop a 
new state preparation routine and show how simple ancilla measurements can efficiently estimate 
the overlap between the linear system solution and an arbitrary state vector. To demonstrate the 
algorithm's applicability, we show how it can be used to compute the electromagnetic scattering 
cross section of an arbitrary target, using the finite element method. This quantum algorithm can 
provide exponential speedup over the best classical algorithm, greatly improving the runtime and 
allowing for the modeling of far more complex objects than possible on a classical computer. 



The potential power of quantum computing was first 
described by Feynman, who showed that the exponential 
growth of the Hilbcrt space of a quantum computer al- 
lows efficient simulations of quantum systems, whereas 
a classical computer would be quickly overwhelmed [1]. 
Shor extended the applicability of quantum computing 
when he developed a quantum factorization algorithm 
that also provides an exponential speedup over the best 
classical algorithm [2]. More recently, Harrow et al. [3] 
demonstrated a quantum algorithm for solving a linear 
system of equations, which in certain cases gives expo- 
nential speedup over the best classical method. 

The importance of linear system solving in scientific 
computing is well known, and the Quantum Linear Sys- 
tems Algorithm (QLSA) represents a potential break- 
through in solving this class of problems using quan- 
tum computing techniques with exponential speedup. 
Prospects for implementation of the QLSA have been 
limited, however, by outstanding issues in state prepara- 
tion (still, in general, an unsolved problem [4-8]), nonuni- 
tarity, and solution read-out not addressed in the original 
paper. These create substantial hurdles in finding appli- 
cations where exponential speedup is possible. 

Here, we present an algorithm that builds on the core 
of the QLSA, but with additions and modifications that 
address the issues of state-preparation, unitarity, and 
read-out necessary to implement the algorithm for prac- 
tical applications. We develop a method to prepare ar- 
bitrary state vectors compatible with the QLSA. Ad- 
ditionally, the QLSA requires a measurement prior to 
accessing information from the solution, preventing the 
ability to use it as a quantum subroutine; we develop 
a technique that avoids use of non-unitary projections. 
To illustrate the functionality of our algorithm, we show 
how one can use it to solve for the scattering cross section 
of an object using the finite element method [9]. Addi- 
tionally, we demonstrate how one can access the solution 
in sub-exponential time, thus preserving the exponential 
speed-up generated by the algorithm. 

In the original QLSA paper of Harrow et al. [3] , the au- 
thors demonstrated a quantum algorithm that could in- 
vert a sparse matrix to solve the quantum linear system 



A\x) = \b), exponentially faster than the best classical 
algorithm in certain situations. The only requirements 
were that the elements of A and \b) can be efficiently 
computable and the matrix A needs to be sparse, or effi- 
ciently decomposable into sparse form. 

Unfortunately, the original algorithm had a few fea- 
tures that made it difficult to apply to a finite element 
scattering computation, and problems of practical inter- 
est in general. Creating the vector \b) generically is an 
unsolved problem, and no mention on how one might 
do this was provided. Additionally their implementa- 
tion was not unitary. The final step required a successful 
outcome of a projective measurement to ensure that the 
quantum state stored the solution. This gave rise to a 
third more subtle issue related to extracting information 
from the solution. Since the solution was stored in a 
quantum state, measurement of it is impractical. For 
most problems of general interest, the QLSA would be 
used as a controlled subroutine, making the original im- 
plementation impossible because of the projective mea- 
surement. Therefore in our algorithm, we remove the 
measurement in order to make the algorithm useable as a 
controlled unitary operator. We anticipate this method 
being useful for other algorithms that make use of the 
QLSA such as the quantum data fitting algorithm [10]. 

The basic outline of the algorithm proceeds as follows. 
First we create the quantum state |6), containing the nor- 
malized right hand side of the linear system, with the 
known values contained in the quantum amplitudes. Al- 
though progress has been made in this area [4-8] , efficient 
quantum state preparation is not in general possible for 
arbitrary states. Within the QLSA we are able to solve 
this by instead preparing the state 

\b T ) = cos&|6)|0)+sin&|&)|l) (1) 

that contains our desired arbitrary state, entangled 
with an ancilla qubit in state |1). This can be done effi- 
ciently in the following manner. We initialize three quan- 
tum registers and an ancilla qubit as 
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An oracle is queried that calculates the amplitude and 
phase components, denoted as bj and <j>j respectively, of 
the vector \b), controlled off the value in the first register. 
Next, apply a controlled phase gate to the ancilla qubit 
in the fourth register, controlled by the calculated value 
of the phase, and finally rotate the fourth ancilla qubit 
conditioned on the calculated value of the amplitude, and 
uncomputc registers 2 and 3 leaving 

I*) -+ "4 E e ^'> (y/^Cffi\0) + C b bj\l)) , (3) 

3=0 

where C'b = l/max(6j) to ensure that all rotations are less 
than 2tt. State (3) is exactly the state (1) with sin 2 fa = 

f EjTo 1 % and cos 2 ^ b = £ ^(1 C 2 b 2 ). This 
technique prepares the desired state \b) entangled to an 
ancilla qubit with value |1) as well as a garbage state, 
denoted by the tilde, entangled to the ancilla with value 
|0). This state can be prepared efficiently so long as the 
unitary operator used to compute bj and (f>j is efficient. 

Next, we apply the QLSA to the state \br)- We modify 
the original algorithm by removing the last post-selection 
step, such that the our implementation is unitary. This 
gives the state 

|*) = (1 - sin 2 fosin 2 ^) 1/2 |4> ) + sin 6 sin X |«) 1 1) 1 1) , 

(4) 

where <f> x a normalization term resulting from the QLSA, 
|$o) is a garbage state in an expanded Hilbert space 
spanned by the solution vector and two ancilla qubits 
which are not in the state |1) simultaneously, and \x) 
is the normalized solution to the linear systems problem 
entangled with two ancilla qubits in the state |1). This 
step is the most computationally expensive part, since it 
requires a Hamiltonian simulation of the matrix A, which 
in its most general form requires 0(k/c) sequential oper- 
ations, where k is the condition number of the matrix A 
and e is the desired accuracy of the solution. The tilde de- 
notes that we have neglected more slowly growing terms, 
similar to the usage in Harrow et al. [3]. 

Next we prepare the state \Rt) = cos<fi r \R)\0) + 
sin</> r |i?)|l) using the same method we used to prepare 
the state \br)- We adjoin this state to Eq. (4) along 
with a fourth ancilla qubit initialized to state |0). We 
apply a Hadamard gate to the fourth ancilla qubit, and 
use it to perform a controlled swap operation between 
the registers containing the solution vector \x) and the 
propagation vector \R), followed by a second Hadamard 
operation on the ancilla. In doing so, we can compute 
the dot product between \x) and \R) as 

IW! 2 = . ^ m °;f m V (5) 
sin (pb sin (p x sin cp r 

where Pi no and Pim refer to the probability of measur- 
ing a 1 in the first three ancilla qubits, and a or 1 in 



the last adjoined ancilla respectively. Because all of the 
techniques used for state preparation and linear system 
solving are unitary, all of the various amplitudes that are 
needed to compute |(P|x)| 2 can be estimated using Am- 
plitude Estimation (AE) [11]. A convenient feature of 
the formalism presented here, is that the Grover oracle 
needed by AE to estimate the various state amplitudes in 
Eq. (5) is trivial to construct. It requires only Z and X 
gates to implement, since the sign change depends only 
on the value of the ancilla qubits (see supplementary ma- 
terial for a detailed explanation). 

We now show how our algorithm can achieve expo- 
nential speedup over the best classical algorithm. On a 
classical computer the runtime is dominated by the linear 
systems solving operation that requires many matrix vec- 
tor products. The best sparse-matrix solving algorithm, 
conjugate gradient, is 0[Ndn log(l/e)] for the scattering 
problem presented here, where d is the number of non- 
zero entries per row and k is the condition number of the 
matrix, while e is the desired precision of the calculation. 

The quantum algorithm requires only O(l) oracle 
queries to create the \b) and \R) states. Estima- 
tion of sin 2 (fit and sin 2 <f> r requires 0(1 /e) iterations 
to estimate to accuracy e with AE. The QLSA re- 
quires Hamiltonian simulation to invert A. Berry et al. 
[12] show show that when using the Suzuki higher or- 
der integrator method [13], this step requires A cxp < 
2m 2 r exp(2-\/ln 5 In (mr/e)) exponential operator appli- 
cations, where m is the number of sub-matrices needed 
to decompose the sparse matrix A into 1-sparse form 
(to = 6gP using the decomposition technique in Rcf. [12], 
where d is the sparsity of A) , and t is the length of time 
the matrix must be simulated times the matrix norm. 
For the algorithm to be accurate to within e, Harrow et 
al. [3] showed that r = 0(n/e). Since we estimate <j> x 
using AE, multiple applications of Hamiltonian simula- 
tion with different times are required. Thus, to estimate 
sin 2 (j) x as well as Pmo and Pun to accuracy e takes 
0(d 4 Klog N/e 2 ) where the tilde again indicates that we 
are neglecting more slowly growing terms. Our imple- 
mentation is quadratically better in k than in the original 
QLSA due to our removal of the post-selection step. 

Combining all steps, the overall quantum algorithm 
has 0(d 4 K,\og N/e 2 ) complexity In situations where the 
condition number scales as poly log(N) the quantum al- 
gorithm is exponentially faster than the best classical al- 
gorithm. However, even if this is not the case, one is not 
restricted to inverting the matrix A as is. Most solvers 
precondition the matrix in some manner to improve the 
conditioning of the matrix. Indeed certain precondition- 
ers are known to reduce the condition number to poly 
log (AT) scaling [14, 15]. In addition, the form of the ma- 
trix A for many practical problems, is known ahead of 
time. This leaves open the possibility for much more 
efficient Hamiltonian simulation than the Trotterization 
method, suggesting greatly improved runtime for specific 



3 



problems, also leading to exponential speedup. 

To demonstrate the algorithm's applicability we now 
show how it can be used to calculate the electromagnetic 
scattering cross section of an arbitrary target. Calcula- 
tion of the scattering cross section is routinely used in 
the electromagnetics modeling community to character- 
ize detectability by radar. In particular, the calculations 
are used to drive design considerations of low-observable 
(stealth) objects. The Finite Element Method (FEM) is a 
numerical technique for solving partial differential equa- 
tions that has been applied to a vast array of problems 
in mathematics, the physical sciences, and other fields 
[9]. In the FEM, a collection of non-overlapping finite 
elements are used to approximate the solution domain. 
These elements are typically simple goemetric shapes 
such as triangles or quadrilaterals in two-dimensional 
spaces and tetrahedra or hexahedra in three-dimensions. 
The fields that satisfy the differential equations are rep- 
resented by a weighted sum of basis functions defined on 
the elements. Along with appropriate boundary condi- 
tions, the problem can then be cast into a linear system 
of equations. 

The FEM can be used as a full-wave solver that con- 
verge to the solution of the continuous problem as the dis- 
cretization length goes to zero. For discretization lengths 
significantly smaller than the length of variation in the 
domain and fields, full-wave solvers provide nearly exact 
solutions. These full-wave methods can be contrasted 
with methods using approximate physics, such as the 
Uniform Theory of Diffraction (UTD) for electromagnet- 
ics [16]. The UTD uses ray optics for propagation and 
scattering from smooth surfaces and accounts for edge 
diffraction by incorporating the analytic solution to scat- 
tering of a plane wave from an infinite wedge. Although 
the UTD provides an excellent approximation for high- 
frequency scattering from simple geometries, it breaks 
down for low frequencies (when ray-optics no longer ap- 
plies) and complex geometries (where simple edge diffrac- 
tion is not the main contributor to scattering). When ap- 
plicable, approximate methods are typically much more 
computationally efficient than full-wave solvers. It is 
clear however, that efficient full- wave solvers are still crit- 
ical for more complex problems. 

The FEM approach to solving an electromagnetic scat- 
tering problem is to break up the computational domain 
into small volume elements and apply boundary condi- 
tions at neighboring elements. This allows one to cast 
the solution of Maxwell's equations into a linear system 

Ax = b. (6) 

The matrix A is constructed from a discretization of 
Maxwell's equation together with appropriate boundary 
conditions due to the scattering object under consider- 
ation. The vector b consists of the known electric field 
components on the scattering boundary. The matrix A, 
which contains information about the scattering object, 



can be efficiently derived from the components of a ma- 
trix F that is dependent only upon the form of the dis- 
cretization chosen to break up the computational domain 
and not the scattering object (see e.g. Ref. [9] and the in- 
cluded supplementary material). Edge basis vectors [17], 
denoted as Nj, are highly popular for electromagnetic 
scattering applications. They give a form of F as 

F t j= f [(V x Ni) ■ (V x Nj) — k 2 Ni ■ Nj] dV (7) 
Jv 

+ ik f (N^-iNj^dS, 

where V is the volume of the computational region, S is 
the outer surface of the computational region, k is the 
electric field wavenumber, the subscript t denotes the 
tangential component, and the indices I and j denote the 
numbering of all the edges contained in the volume V. 
The surface integral is an absorbing term used to prevent 
reflections off the artificial computational boundary. For 
simplicity we assume a first order boundary condition; 
however, that limitation is not required for the quantum 
algorithm. On the inner scattering surface the correct 
boundary condition for the scattered field on metallic 
scatterers is n x E = — n x EW, where EW is the in- 
cident field, E is the scattered field, and n is the unit 
vector normal to the surface is applied. The edge basis 
vectors are defined such that V • Ni = 0, which implies 
V • E = and constant curl, making them valid basis 
functions to use for the electric field. A much more de- 
tailed explanation of the FEM construction can be found 
in the supplementary material. 

The actual electric field at any point can be con- 
structed from the linear system solution by interpo- 
lating with the edge basis vectors giving E^ e '(r) = 
N^x\ e \ Here the superscript (e) denotes the par- 
ticular element in which the field is being calculated, and 
n is the number of edges of the element. The radiated 
electric field in the far-field at position s is given by 




+ iks xM(r)} e ikSr dS, 



where J(r) = — (iw/io) _1 n(r) x V x E(r) and M(r) = 
E(r) x n(r), with n(r) the unit normal of the surface S 
at r. Using the edge basis expansion gives 

E(s)-p= f — J2Rk(s)x k , (9) 

k 

where p is the radar polarization (with p • s = 0) and 

i? fe (s,p) = p- / s x {s x [(V x N fc ) x n] (10) 
Js 

+ ikN k x n}e mr dS, 
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where the index k here is the global edge index. The 
radar scattering cross-section (RCS) in the direction s is 
given by 



RCS = lim 47rs 2 |E(s) • p| 2 = -*-|R • 



(ii) 



or simply the dot product of R with the solution x, where 
we have assumed an incident plane wave with unit electric 
held amplitude without loss of generality. 

The edge basis elements can take a simple functional 
form, which allows one to analytically evaluate the in- 
tegrals in Eqs. (7) and (10). This allows for efficient 
computation of the matrix and vector elements, a re- 
quirement for the quantum algorithm. Because of the 
local nature of the finite element expansion, the volume 
and surface integrals extend only over the region encom- 
passed by the finite element. As a result A is highly 
sparse, allowing an efficient decomposition into a 1-sparse 
form [12], also necessary for the quantum algorithm. 

To obtain the cross section using the quantum algo- 
rithm one uses the oracles just presented to create the 
A matrix and |&) and \R) state vectors. Then one must 
restore units to the normalized output received from the 
quantum algorithm. Doing so yields the following equa- 
tion for the cross section in terms of outputs from the 
quantum computation 



RCS = 



1 N 2 



t>b sin 



47T 



C 2 C 2 sin 2 ■ 



-(Pmo -Pun), (12) 



where Cb — l/max(b) and C r = l/max(R) are known 
parameters. Thus to compute the cross section, we esti- 
mate each sin 2 4>(b,x,r) term as well as the Pmo and Pun 
terms independently using AE. 

Finally we remark on the efficiency of the scattering 
cross section calculation. With no preconditioning, fi- 
nite element condition numbers scale as N 2 / n [18, 19], 
where n is the number of dimensions of the problem, im- 
plying that even in the most general case our algorithm 
scales better than its classical counterpart for a three- 
dimensional finite element problems. However, as men- 
tioned above there are likely quantum preconditioners or 
better simulation techniques that can further improve the 
quantum scaling to provide exponential speedup. 

We have demonstrated a quantum algorithm that gen- 
eralizes the QLSA to solve an arbitrary linear systems 
and estimate the overlap of the solution with an arbi- 
trary state vector by simple ancilla amplitude measure- 
ments. To demonstrate its functionality we showed how 
one could use it to solve an electromagnetic scattering 
problem using the finite element method and estimate 
the scattering cross section. In the best case, where the 
Hamiltonian simulation or condition number scaling can 
be reduced to poly log (A), this algorithm has complexity 



that is exponentially better than the best classical algo- 
rithm. This opens up the potential for quantum comput- 
ing to be applied to a broad class of problems of practical 
interest to the computational physics community. 
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